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We derive and employ a semi-classical Langevin equation obtained from path-integrals to describe 
the ionic dynamics of a molecular junction in the presence of electrical current. The electronic envi- 
ronment serves as an effective non-equilibrium bath. The bath results in random forces describing 
Joule heating, current-induced forces including the non-conservative wind force, dissipative frictional 
forces, and an effective Lorentz-like force due to the Berry phase of the non-equilibrium electrons. 
Using a generic two-level molecular model, we highlight the importance of both current-induced 
forces and Joule heating for the stability of the system. We compare the impact of the different 
forces, and the wide-band approximation for the electronic structure on our result. We examine 
the current-induced instabilities (excitation of runaway "waterwheel" modes) and investigate the 
signature of these in the Raman signals. 



PACS numbers: 85.75.-d, 85.65. +h,75.75.+a,73.63.Fg 

I. INTRODUCTION 

The interaction of electrons with local vibrations 
(phonons) has an important impact on the conduction 
properties and stability of molecular conductord^^ anc j 
has undergone intense study both experimentally and 
theoretically^^. In the low bias regime where the volt- 
age is comparable to phonon excitation energies, valu- 
able information about the molecular conductor can be 
deduced from the signature of electron-phonon interac- 
tion, known as inelastic electron tunneling spe ctroscopy 
(IETS), or point contact spectroscopy (PCS) 2 ' 15 ' 49 ^ . 
Theoretically, this regime has been addressed with some 
success u sing me an- field theory such as density functional 
theory^! ! 23 * 37 ! 49 !, where the vibrations are assumed to 
be uncoupled to the electrons, while the effect of the 
phonons on the electronic transport is taken into ac- 
count using perturbation theory. On the other hand, for 
higher voltage bias, and for highly transmitting systems, 
a large electronic current may strongly influence the be- 
havior of the phonons even for relatively weak electron- 
phonon coupling. The resultant "Joule he ating" is well- 
known in the molecular electronics contexlP 11 ! 12 * 26 !, and 
remains a lively area with a range of approaches (and 
with occa siona l lack of complete agreement between 
treatments" 251 ). 

More recently the current-induced wind-force known 
from electromigratioiP 2 ' has been reexamined for atomic- 
scale conductors and shown also to be able to excite the 
conductor and possibly lead to a runaway instability^^!. 
It has been shown how a part of the force on an atom in 
the presence of the current may have a non-conservative 
(NC) component, able to do net work around closed 



paths. This was explicitly proven by calculating th e curl 
of the vector field describing the force on an atorrPSMl. 
The NC energy transfer - also dubbed the atomic "wa- 
terwheel" effect - requires a generalized circular motion 
of the atoms and involves the coupling of the electronic 
current to more than one vibrational mode. 

Along with the NC force contribution we have re- 
cently identified a velocity-dependent current-induced 
force which conserves energy and acts as a Lorentz-like 
force on the generalized circular motion. This force can 
be traced back to the quantum mechanical Berry-phase 
(BP) of the electrons^ES. Together with the NC force we 
will, further, have a component which is curl-free and is 
related to the change in the effectiv e pote ntial energy sur- 
face of the atoms due to the currenlP^J as well as nona- 
diabatic "electronic friction" forces^. In the nonequilib- 
rium situation the "friction" force can, however, turn into 
a driving force amplifying the vibration. This happens 
under certain resonance conditions akin to a l aser e ffect, 
but now involving phonons instead of photon^EH 

A unified approach including all aforementioned effects 
on an equal footing is highly desirable for further study 
in this direction. In this paper, we extend the electronic 
friction approach proposed by Head-Gordon and Tully^l 
for molecular dynamics to take into acco unt the no nequi- 
librium nature of the electronic currenlP^ 59 1 61 1 641 ^. A 
similar approach has been taken to descr ibe m odels of 
nano-electro-mechanical systems (NEMS| 67 * 68 l Using 
the Feynman- Vernon influence functional approach, we 
derive a semi-classical Langevin equation for the ions, 
which we can use to study Joule heating, current-induced 
forces, and heat transport in molecular conductors. We 
perform a perturbation expansion of the electron effec- 
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tive action over the electron-phonon interaction matrix. 
This allows us to make connections with other theoretical 
approaches, especially the nonequilibrium Green's func- 
tion (NEGF) method, used to study the Joule heating 
problem* 20 * 23 *. We also give an extension of the perturba- 
tion result to the adiabatic limit, which makes connec- 
tions with our previous results, and solves an infrared 
divergence problem in the expression of the BP force in 
Ref. [55] We apply the theory to a two-level model in 
order to (1) clarify the roles played by different forces re- 
garding the stability of the device, and (2) discuss the sig- 
nature of the current-induced excitation in the Raman- 
scattering especially focussing on conditions close to a 
current-induced runaway instability. 

The paper is organized as follows. In Sec. [IT] we briefly 
review the derivation of the generalized Langevin equa- 
tion. In Sec. |III|we analyse the electronic forces entering 



into the Langevin equation. Sec. IV compares the effect 
of different current-induced forces for a two-level model, 
concentrating on the NC and BP force. In Sec.|V| we ex- 
tend the perturbation result to the adiabatic limit, and 
introduce coupling of the system with electrode phonons. 
The derived formulas can be used to study the current- 
induced phononic heat transport. In Sec. |VI[ we present 
ways of calculating the quantum displacement correla- 
tions, which is essential for the theoretical description 
of Raman spectroscopy in the presence of current. Sec- 
tion |Vn] gives concluding remarks. 



II. THEORY 
A. Influence functional theory 

We start from the influence functional theory of Feyn- 
man and VernorP^, which treats the dynamics of a "sys- 
tem" in contact with a "bath" or "reservoir" . In our 
case the system of interest consists of the few degrees 
of freedom describing the ions of a molecular conduc- 
tor, interacting with the electronic reservoir composed 
of all electronic degrees of freedom in the molecule and 
electrodes, as well as the phonon reservoirs of two elec- 
trodes. The reservoirs can be out of equilibrium gen- 
erating an electronic current, and may, further, involve 
a temperature difference generating a heat flux between 
the electrodes. All effects of the bath are included in the 
so-called influence functional, which gives an additional 
effective action, modifying that of the isolated system. 

Now we briefly review the idea of the influence func- 
tional approach. With the help of the influence func- 
tional, F, the reduced density matrix of the system in 
the displacement representation reads, 



<^2|p s (*2)|2/2> = / dxi J dyifC(x 2 ,y2;xi,yi) 

x(xi\ps(h)\yi), (1) 



with the propagator of the reduced density matrix 

l-{x,v)(t 2 ) = (x2,V2) 



V{x,y) 



x exp 



-(s s (x)-s s (y)) 



Here x and y are a pair of displacement histories of the 
ions, and S s is the action of the system only. In deriv- 
ing this, we have assumed that the system and bath are 
uncorrelated at tt(—> — oo), 



p(ti) = p a (h) <g> Pb(h). 



(3) 



The influence functional includes the information of the 
bath S b and its interaction with the system S 1 ,, 

r r(r,q)(t 2 ) = (r 2 ,r 2 ) 

F(x,y)= dr 2 dndqi / V(r,q) 
J •'(r,9)(ii)=(n,5i) 



xexp 



(S b {r) + Si{x, r) - S b (q) ~ Si(y, q)) 



x{ri\pb(ti)\qi) 



(4) 



with r and q representing forward and backward paths of 
the bath degrees of freedom. Most importantly, a correc- 
tion to the action of the system can be defined from the 
influence functional AS = —ih hiF(x,y), which is usu- 
ally not time-local or real. It has been used to derive a 
semiclassical Langevin equation, describing the dynam- 
ics of the system interacting with the environment!^^. 
In this case new variables are introduced, 



Q = ^{x + y), Z = x-y, 



(5) 



describing the average and difference of the two paths, 
respectively. In the semi-classical approach the average 
path, Q, is shown to yield the variable in the Langevin 
equation, whereas role of the difference, £, is to introduce 
fluctuating random forces in a statistical interpretation. 
The influence of the environment will favor paths with 
small excursions given by £, and will ensure that only the 
solution, Q, obeying the classical path will contribute for 
a high temperature reservoir. We illustrate this further 
below. 

Next, we introduce our model and give the result for 
the influence functional describing the nonequilibrium 
electron bath. From the effective action, we can read 
out the forces acting on the ions due to the electrons. 
We also discuss how a thermal flux may be included. 
Parts of th e deriva tions can be found in our previous 
publication d 55 * 64 *^, but here we aim at a more general 
formulation, which we present with detailed derivations 
together with an illustrative model calculation. However 
we note that the theory is fully compatible with more 
realistic systems with complex electronic and vibrational 
structure treated within a mean-field approach such as 
density functional theory. 
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B. System setup and Hamiltonian 

To obtain an effective action describing the vibrations 
(phonons) in the molecular conductor we first divide 
the complete system into electron and phonon subsys- 
tems. We will treat electron-electron interactions at the 
mean-field level. To describe the nonequilibrium situa- 
tion where a current is flowing through the molecular con- 
ductor between two reservoirs, the electron subsystem is 
further divided into a cental part (C) and two electrodes 
(L,i?), whose electrochemical potentials change with ap- 
plied bias. For the purposes of the present study, we allow 
electrons to interact with phonons in C only, and further- 
more ignore the anharmonic coupling between these dif- 
ferent modes. The coupling of the molecular vibrations 
with electrode phonons will be considered in Sec. |VB| 

The single particle mean-field electronic Hamiltonian 
at the relaxed ionic positions, H , is written within a 
tight-binding or LCAO type basis with corresponding 
electron creation (annihilation) operators for the jth or- 
bital, Cj (cj). The Hamiltonian H° spans the whole LCR 
system, while the electron-phonon interaction is localized 
in 6^22]. The Hamiltonian of the whole system reads, 



H = 



H. 



ph 



Hph 

1 rp 
2 P P - 



H e {u), 
1 T 

—u Ku, 
2 



(6a) 
(6b) 



H e {u) = 



ij C i C j 



E 

k,i,j€C 



M*c\ Cj u k , (6c) 



where p are momenta conjugate to u, and u is a col- 
umn vector containing the mass-normalized displacement 
operators of all ionic degrees of freedom (e.g. Uk = 
^mk(rk ~ r k)i where mk is the mass of ionic degree of 
freedom k and (rj?) is its (equilibrium) position). The 
equilibrium zero-current dynamical matrix is denoted by 
K. This Hamiltonian has been used to describe IETS in 
molecular contacts with param eters obtained e.g. with 
DFT for concrete systems^^l. We have previously dis- 
cussed the adiabatic limit, where the perturbation is in 
terms of the velocities ilk and where the full non-linear 
effects of Uk can be included. Here we instead assume 
small displacements from equilibrium and expand the 
electronic Hamiltonian to first order in Uk- Later in 
Scc. |V A we compare this to the adiabatic limit, discussed 
in Rets. 1551591 

The correction to the action due to the coupling to the 
electron reservoirs can be found using the linked-cluster 
expansion in the coupling, M fc , following Ref.1551 The ef- 
fective action of the nonequilibrium, noninteracting elec- 
tron bath reads, 



AS{x,y) = ihY] f dX f dr 
■ Jo Jk 



x Tr[G(T,T + ;X)M k X k (T) 



(7) 



where the trace Tr [] is over the electronic bath and this 
will be so in all following formulas. The parameter A 



is used to keep track of the order in the linked-cluster 
expansion^. The time r is defined on the Keldysh 
contour--' K. On the real time axis the Green's func- 
tion decomposes into 



G(r,r') = 



G(t,t') G<(t,t') 
G>(t,t') G(t,t') 



and 



X(T) = ( *(*> 





-y(t) 



(8) 



(9) 



Time r+ is infinitesimally later than r on the whole 
Keldysh contour. The limits of integration extend to -co 
and +oo if not specified. This applies to all the integrals 
in the paper. The Green's function is given by the Dyson 
equation, 

G(t,t + ) = G q (t,t + ) (10) 

+ E / dT'G Q (T,T')M k X k (T')G(T>,T + ), 
k JK 

Gq being the single electron Green's function without 
interaction with the ions, which reads 



(11) 



2ttH 

np(e — n a ) — 9(t — t 1 ) np(£ — /i Q ) 

np(e-fi a )-l n F (e - (j, a ) - 6{t' - t) 



n F {e - n a ) = 1/ [l + e ( E -^)As T ] is the Fermi-Dirac 
distribution function for electrode a, A — ^ Q A a , and 9 
is the Heaviside step function. The spectral function is 
defined as 



^(e)=iG5( e )[E;( e )-E«( e )]Gg( e ), 



(12) 



and S^(e) (E^(e)) the retarded (advanced) electron self- 
energy from lead a. We use e and w as parameters for 
the electron and electron-hole pair/phonon properties, 
respectively. 

The effective action in Eq. (7| could be expanded into 
an infinite series in M. We only keep terms up to second 
order, assuming small M (or alternatively assuming small 
displacements as stipulated earlier). The first-order con- 
tribution written in terms of the average and difference 
paths for the vibrations, Qk,£,ki reads, 



(13) 



with a first-order, displacement-independent force term, 



-2E 



de 
2n 



TiiA a (e)M k ]An F (e), (14) 



An F (e) = n F (e - jj, a ) - n F (e - /U ) 



(15) 



Here (Iq is the equilibrium electrochemical potential and 
the factor of 2 accounts for spin degeneracy. Above we 
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explicitly subtract the equilibrium forces (obtained for 
Ma = Mo), since these are already included in the elastic 
forces described by K. The different filling of electronic 
states originating from different electrodes in Eq. (fl5| re- 
sults in a displacement-independent "wind" forc cP^ 1 '^ 1 ^ 1 . 
Its effect amounts in the small-displacement approxima- 
tion used here to a bias-induced shift of the equilibrium 
ionic positions. We will therefore ignore it from now on, 
since we will be considering the effects of the nonequilib- 
rium electrons on the ionic dynamics. 

Before introducing the second-order contribution, we 
note that the applied bias between the two electrodes 
also modifies the electronic Hamiltonian, and thus A a 
and M k . It is the result of charge rearrangement in the 
device in response to the applied bias. This is out of the 



scope of present paper, and is not included in Eq. (15). 
The inclusion of external electric fields in the electronic 
Hamiltonian can account, for example, for the "direct" 
electromigration forced. 

We now turn to the second-order contribution central 
to our discussion, 

AS (2) (<3,£) = -~ J dtj dt' J du 



x e -*"(*-*') A°f(w) 



a,p,l,k ' 



coth 



hw - (n a - flp) 



2k B T 



26{t' -tMt')Q k {t) 



(16) 



where the central quantity is an interaction-weighted 
electron-hole pair density of states (incl. spin), 



x Tr[M k A a (e 1 )M l A p (£ 2 )] 

x (n F (ei - fi a ) -np(e 2 - Hp)). 

It has the following properties, 

A2f( W )=Agf 

A2f( W ) = -A^(- W ) ) 
A(-w) = -A», 

where we have defined A = ^ Q « A"^. 



(17) 



(18) 
(19) 
(20) 



C. The generalised Langevin equation 

In order not to complicate the equations we will in the 
following suppress the phonon-mode index and implicitly 
write vectors and matrices without these. Note that these 
phonon indices are generally not interchangeable. This 
calls for care for example when carrying out permutations 
within the trace in A?, and quantities derived from it. If 



the reader wishes to make such rearrangements, it is nec- 
essary to reinstate the indices (fc, / above), in the correct 
starting order, first. Using the results of Subsec. |II A| wc 
get the path-integral, 



K 



Vi \ VQexp 



dt / dt'£ 2 (t)L(t,t')Q(t') 



x exp 



1 

2h 



dt / dt'C (t)u(t - 



(21) 



We have defined 



and 



with 



IT(i-t') = 2mJ- 1 {A(uj)}, 
W(t-t') = 6(t-t')I[(t-t'), 



n(t-f)=j- 1 {nM}, 



(23) 
(24) 



(25) 



n(w) = HoH + An(u) 

= — 7rA(w)coth ( z—= 



-TT^A^M 

coth ( ^-^-^ ] ) coth f hL ° 



2k B T 



\2k B T 



(26) 



VB 



We have split II into two terms. We will see in Sec. 
that All is responsible for the heating effect. The 
Fourier transform is defined as J-{f(t)} = J dtf(t)e tu:t , 
and J-H/M} = / f^/Me-^. After a Hubbard- 
Stratonovich transformation, we get 



K 



Vi\VQ\ Vf 



x exp 

x exp 
x exp 



^)exp(-^i 



dt^it) / dt'L(t,t')Q(t')-f(t) 



2h 



dt dt'fm-^t^fif) 



(27) 



The factors in line 2 (where £i = £(t%), etc) arise from 
the integration by parts to transform the kinetic energy 
(£ T (t)Q(t)) into the form in line 3, and are needed for 
example if one wishes to make a connection with the 
Wigner function. The above form of the effective action 
suggests a classical interpretation to the motion of the av- 
erage di splace ment Q. It follows a generalized Langevin 
equatiorPEI] 

Q{t) = -KQ(t)- flL r (t-t')Q(t')dt' + f(t), (28) 
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where f(t) is a classical stochastic force, whose time- 
correlation is given by 



properties: 



(f(t)f T (t')) = nii(t-t'), 



(29) 



Equation ( 28 1 is the equation of motion for harmonic os- 



cillators, perturbed by the second and third terms due to 
interaction with electrons. There are two ways of seeing 
how it arises. First, it is the Euler-Lagrange equation of 
motion for the action in Eq. ( 27 ). Alternatively, if in ( 27 ) 



we carry out the ^-integral, that introduces 5(L ■ Q — /) 
(where • denotes time-integration for short) in the inte- 
grand, which restricts Q(t) to the evolutionary path gen- 
erated by (28 1. The physical significance of Q(t) is that it 
is the quasi-classical coordinate appearing in the Wigner 
function. Its Newtonian equation of motion above relies 
on the quadratic nature of the effective action as a func- 
tional of the generalised coordinates, and would have to 
be revisited in the presence of higher-order terms in the 
Hamiltonian, or in the Green's function expansion. 

It is possible to solve the generalised Langevin equa- 
tion by Fourier transform. From that, we get the semi- 
classical displacement correlation function, 

\ (QQ T )(u) = D r (uj)tl(uj)D a (uj). (30) 
n 

Here the phonon retarded Green's function is defined as 

D r (uj) = {D a (Lu))l = [(u + i0 + ) 2 - K - IT»] -1 . 

(31) 

Note that II and II can be written as the standard phonon 
self-energies in the NEGF method 

n(w) = n r (u;) - n a H = n>(w) - n<(w), (32) 



fiH 



[n>H + n<H] 



(33) 



The self-energy diagram is shown in Fig. [TJ and their 
expressions are given in Appendix [Cj 




FIG. 1: Lowest order phonon self-energy due to interaction 
with electrons. 



III. FORCES 

A. General results 

The electronic forces in the Langevin equation are di- 
vided into stochastic and deterministic parts. The corre- 
lation function of the stochastic force has the following 



nt(«) = n(w), 
rt(-w) = ft*(w). 



(34a) 
(34b) 



Consequently, U(t) is real, but in general fi(— t) ^ Yl(t) 
at finite bias. 

Now let us look at the deterministic forces due to elec- 
trons. TI(uj) has the properties: 



nt(«) = -fi(w), 
fi(-w) = n*(cj). 



(35a) 
(35b) 



In the Fourier domain, we can split it into different con- 
tributions, 

-lT(u;) = -i7rReA(w) + 7rImA(w) 
-7r-H{ReA(w')}(w) - wr^{ImA(u/)}(u>), (36) 

where the Hilbert transform is defined as H{g(x')}(x) = 
—V j tril dx'. Using the symmetry properties of A(lo), 



we can now examine each term in Eq. ( 36 1 



The first term is imaginary and symmetric. It de- 
scribes the standard friction, i.e. processes whereby the 
motion of vibrating ions generate electron-hole pairs in 
the electronic environment. This process exists also in 
equilibrium. We can write 



Ffr = -T)(u))Q(u)), 
with the friction matrix defined by, 



r)(u>) 



-ReA(w). 



(37) 



(38) 



The second term in Eq. ( 36 1 is real and anti-symmetric. 



It has a finite value even in the limit of zero fre quency. 
It is describing the NC force, discussed recently^ 156 1 58 1 59 1 



with 



F NC =Af(u)Q(u), 



Af(uj) = 7rImA(cj). 



(39) 



(40) 



The third term is real and symmetric, and can be con- 
sidered a renormalization (RN) of the dynamical matrix 



with 



f rn = -CMQH, 



CM = nH{ReA(oj')}(ui). 



(41) 



(42) 



Finally, the last term is imaginary and anti-symmetric, 
proportional to u) for small frequencies. Hence it is to be 
identified with the BP force in Ref. [551 



F BP = -B(u)Q(u), 
with the effective magnetic field 



B(u) = H{ImA(u')}(u). 

CO 



(43) 



(44) 



G 



B. Equilibrium and nonequilibrium contributions 

We can divide the A(w) into an equilibrium part and 
a nonequilibrium part A(lo) — K eq {uj) + A A (a;), and look 
at their contribution to the forces separately. A eq (uj) is 
given by Eq. ( fl7| ) with /i Q and [ip replaced by the equilib- 
rium electrochemical potential no- The nonequilibrium 
part can be written as 



Im 

Re 



)aAM = 2£ f A A »«( e ) (45) 

a J 

MA a (e)M ( A(e.) ( + )a(e+)] 



Im 

Re 



Tr 



with e± = e ± fiw. 

In the following we assume zero magnetic fields and 
treat M k and A(e) as real symmetric matrices in the 
electronic real-space basis. 



1. Equilibrium contribution 

We consider first the equilibrium part A eq (u>). It is 
real, giving the equilibrium friction, and its Hilbert trans- 
form gives the equilibrium renormalization of the poten- 
tial. 

Friction - The equilibrium friction matrix reads 



1 



de 



n F {e 



Mo J 



eqK ' 2oj J 2tt' 
x Tr [MA(e)M (A(e+) - A(e-))] . 

Renormalization - The equilibrium RN reads 
de 



Ce 9 (w) = 2 



2tt 



n F (e - no) 



(46) 



(47) 



x Tr [MA(e)M (R(e-) + R(s+))} 
We have defined 



R(e) 



l -U{A{e')}{e) 



G r (e) 



G Q (£) 



(48) 



In general Ceg(w) has a frequency dependence, of 0(uj 2 ) 
or higher. Its static (frequency-independent) part is al- 
ready included in the dynamical matrix, when calculated 
within the Born-Oppenheimer approximation. 



2. Nonequilibrium contribution 

Now we consider contributions from the the nonequi- 
librium part AA(w). Although we are considering the 
two-terminal LCR assembly described earlier, an arbi- 
trary number of independent terminals can be accom- 
modated via the summations over indices a, /?. In the 
two-terminal case, we write eV = fj,L ~ Mfl where V is 
the bias. 



Friction -We first get a correction to the equilibrium 
friction 



a " 

x ReTr [M(A(e. 



de An F (e) 



2lj 

A{e-))MA a {e)\ 



(49) 



This nonequilibrium correction may give rise to another 
interesting instability, characterized by a negative fric- 
tion, if the spectral functions depend on the energy 
in a special way which enable a population-inverted 
situational. It is responsible also for enhanced heating, 
and for the converse: current-facilitated thermal relax- 
ation, in systems with appropriate spectral featured 76 ! 77 !. 
NC force - The coefficient for the NC force is 



Af(u,) = 2£ J 



de An£,(e) 



2tt 2 

x ImTr [MA a {e)M{A{e+) + A(e_))] 



(50) 



Performing the Hilbert transform, we get the nonequilib- 
rium correction to the RN force and the BP force. 

Renormalization - The nonequilibrium correction to 
the RN force is given by the coefficient 

ach = 2]T / Ja<( £ ) 

a 

x ReTr [MA a (e)M(R(e + ) + R(e-))] . (51) 



BP force - The BP force is 



F bp (lo) = -B(w)Q(w), 



(52) 



with 



de An F (e) 
2tt lo 



x JmTr[MA a (e)M(R(e+)-R(e-))]. (53) 

The relative magnitude of the BP force and NC force can 
be estimated as, 



IF, 



BP 



I -Five | 



ftuj 

W\ 



(54) 



\H\ being a typical electronic hopping integral. Ref. 1551 
instead of \H\, used the phonon frequency as a cutoff 
when calculating the BP force. This over-estimates the 
effect of the BP force. A more detailed discussion of this 
is given in Appendix \K\ 

If the dynamics of the ions is very slow compared to the 
dynamics of the electrons, the electronic spectrum varies 
weakly within the vibrational energy spectrum, and we 
can take the U) — > limit in the expressions for the forces. 
All deterministic forces then become time-local. This will 
be compared with the adiabatic result in Sec. V A 
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C. Wideband approximation 



2. Phonon excitation 



A regime of practical interest is the limit where the 
electronic spectrum varies slowly not only within the vi- 
brational energy spectrum, but also within the bias win- 
dow. Then we can ignore its energy dependence alto- 
gether and evaluate all electronic properties at the Fermi 
level (no)- This is the wideband approximation used for 
example in Ref. OH to study IETS. 



1. Dynamical equations 

Within the above approximation, the Langevin equa- 
tion reads 



Q(t) = -KQ{t)- m Q(t)+M Q{t) 
- ( Q(t)-B a Q(t) + f(t), 



with 



Vo 
•A/o 

Co 



2— Tr[MA( Mo )AM( M o)], 
eV X -, 

2^ReTr [MAA((i )MR([i )] 
2ir 



o 



heV 
'~2m 



ImTr [MAA(iJ,o)Md s R(no)} 



where we have introduced, 

x - = 2^ImTr[MA L (no)MA R (v )], 
X+ = 2±-ReTr[MA L (» )MA R ( tM) )], 



(55) 

(56) 
(57) 
(58) 

(59) 

(60) 
(61) 



and Ayl(/io) = Al(/.iq) — Ani^o)- The noise correlation 
function also takes a simpler form, 



n (w) = (ujt] - ieVx ) coth 



\2k B T 



and 



An(w) 



(62) 



(63) 



coth 



/ huj + aeV 
\ 2k B T 



coth 



\2k B T 



The noise originates from the fluctuating part of the 
forces (Eq. ( |B4[ )), whose correlation spectrum in fre- 
quency space is in general Hermitian, but not real at 
finite bias. The real part corresponds to the friction, and 
the imaginary part to the NC and BP force in the de- 
terministic forces. Importantly, the quantum zero-point 
fluctuation is taken into account (encoded in the coth 
function). 



When isolated from the electrode phonons, the system 
could get heated or cooled due to the passing electrical 
current. At steady state, the system phonon population 
will be different from that at equilibrium. We will now 
compare the phonon excitation result from the Langevin 
equation with that from NEGF theory^ 1 . To do that, 
we employ the wideband approximation, and ignore cou- 
plings between different phonon modes. In this way we 
can study each mode separately. At steady state, the en- 
ergy stored in each phonon mode can be obtained from 
the solution for Eq. ( 55 ) in frequency space, 



Ei = (<#(*)) 

= / u?(QiQi)(w) 



2tt' 



(64) 



Using Eq. ^ and Eqs. ( |62p3| ) , assuming small broad- 
ening of the phonon mode, it can be further simplified 

as 



~2~ C ° \2k B T 



= IN, 



1 



(65) 



Here we have introduced an effective phonon number N%. 
At low temperature AIIjj(a;) can be approximated as, 

AUuiuj > 0) oc (eV - hw)0(eV - Huj). (66) 

Figure [i] shows the bias dependence of AII^ at zero and 
finite temperature for a phonon mode huji =0.1 eV. 
Interestingly, the Joule heating exhibits a threshold for 
phonon excitation at the phonon energy, at zero tem- 
perature. If fact, Eq. (65) is exactly the same as the 



quantum result Eq. (47) in Ref. 20. If we take the dif- 
ferent definition of the electron-phonon interaction ma- 
trix used here and in Ref. HOI we find that the friction 
coefficient rju, and the nonequilibrium noise spectrum 
AILii(oJi) are related to the electron- hole pair damping 
Jl_ h , and phonon emission rate 7* TO defined in Ref. [2131 as 

ll-h = Vu, AftiifWi) = 2ucf\ m . We therefore conclude 
that we recover the quantum-mechanical result from the 
semi-classical Langevin equation. This is because we 
only need the quantum average of the equal-time dis- 
placements (ui(t)uj(t)) to study the energy transport 
which can be calculated exactly from the semi-classical 
Langevin equation (see Appendix [B] for details). Alter- 
natively, in Wigner-function language, the mean phonon 
energy is expressible solely in terms of the coordinate Q 
(and/or the velocity Q). As we have seen, for a harmonic 
action the present Newtonian equation of motion for Q 
is exact. 




is defined as eV = \il — and the average electrochemical 
potential fj, = (fiL + Mfl)/ 2 - 



FIG. 2: Example of the nonequilibrium noise spectrum at 
T = K (solid) and 300 K (dashed) for a given phonon mode 
with frequency Hljo = 0.1 eV. 



IV. NUMERICAL RESULTS FOR A SIMPLE 
TWO-LEVEL MODEL 

To build an intuitive understanding of the theory 
above, we now apply it to a simple spinless two-level 
model which could describe a diatomic molecule. For 
this model system, it is possible to do the calculation 
using the general results in Sec. |III A| without the ap- 



proximations developed in Sec. |III C We start from a 



model Hamiltonian for the isolated system, 

H = H e + H p h + H int 

= e (c|ci - c\c 2 ) - t(c\c 2 + cjci) 



(67) 



The electrons couple with two phonon modes in the fol- 
lowing two forms: 



and 



H Lt = miui(c\c2 + cjci), 



H int = m 2 u 2 (c\c 1 - c\c 2 ). 



(68) 



(69) 



The first and second electronic levels couple with the left 
and right electrode, respectively, with level broadening 
r. A sketch of the model is shown in Fig. [3j Mode 1 
corresponds to the bond-stretching mode of the diatomic 
molecule, while mode 2 mimics the rigid motion of the 
diatomic molecule between the two electrodes. 



A. Current- induced forces and phonon excitation 

Figure[4]shows different parts of A(w) and their Hilbert 
transforms. The solid and dashed lines in the bottom- 
right panel correspond to the NC and BP force. For the 
parameters used here, the NC and BP force are compara- 
ble with the diagonal RN and friction forces, shown in the 



two panels on top. The symmetry properties imply that 
the NC and the RN term become dominant in the limit 
of slow vibrations. The nonequilibrium RN term can be 
of interest , for example qualitatively changing the poten- 
tial profile^U£7'. Furthermore, by combining the RN term 
with a further contribution, arising from the next order 
in the expansion of the electronic Hamiltonian in powers 
of the displacements, it is possible to construct the full 
non-equilibrium dynamical response matrix^. Its bias- 
dependence can compete with the NC force and influence 
the appearance, or otherwise, of waterwheel modes. 

In the present work the RN term only changes quan- 
titatively the results for the model used. The RN term 
will be excluded altogether below, focusing instead on 
the effect of NC and BP force. 

We already see from the analytical result that the mag- 
nitude of the BP force is directly related to the energy 
dependence of the electron spectrum. This is confirmed 
numerically in Figs. [5] and [6] We show in Fig. [5] the rela- 
tive magnitude of the BP force compared with the aver- 
age diagonal friction for different level broadenings T at 
1 V with hu>o — 0.02 eV. The inset shows the left spectral 
function. We see that for a range of T, the BP force is of 
the same magnitude as the friction. With increasing T 
the resonance in the DOS gets broader, and consequently 
the BP force gets smaller (oc d £ R as in Eq. ([59}). In Fig. D 



we vary the energy position of the bias window relative 
to the peak in the spectral function. The BP force drops 
quickly when the bias window moves away from the peak. 

Assuming a small detuning of the two harmonic oscil- 
lators hw± = hujQ ±6/2, we now study the bias depen- 
dence of their frequency, and damping described by their 
inverse Q-factors in Figs. [7] and [8] The runaway solution 
is defined at the point where the damping disappears, 
1/Q» = — 2Ima;i/Rea;i = 0. We see that the BP force 
in general helps the runaway solution by reducing the 
threshold bias. This is prominent for larger detuning. 
The reason is that it bends the eigenmodes into ellipses 
so that the NC force continuously can take energy out of 
one mode, while pumping energy into the other. Eventu- 
ally, this changes the polarization of the harmonic motion 
from linear to elliptical (circular) in mode space. 

Comparing the full calculation with that from the 
wideband approximation, we see that they agree well 
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FIG. 4: Different parts of the A(u>) function (solid) and their 
Hilbert transform (dashed). The model parameters are F = 1 
eV, t — 0.2 eV, Eo = 0, mi = = 0.01 eV/y / amuA, /io = 0, 
and V = 1 V. Indices 1 and 2 refer to the respective phonon 
modes. The equilibrium renormalization term has already 
been substracted in the plot. We use the same parameters in 
the following figures if not stated explicitly. 



only in the low bias regime. For large bias, we need 
to take into account the energy dependence of the elec- 
tronic spectral function. The frequency dependence of 
the threshold bias with and without the BP force is de- 
picted in Fig. [9] The divergent behavior of the threshold 
bias when only the NC force is considered is due to the fi- 
nite range of the electron DOS (inset of Fig. [5]) . Once the 
bias is large enough for the bias window to enclose the 
DOS peaks, the NC force will saturate. Further increase 
of the bias does not help. But the BP force has an extra 
linear u> dependence (since it depends on Q, instead of 
Q), which becomes important for high frequencies. The 
bias dependence of the mode-correlation function and de- 
rived excited phonon number (cf. Eq. (64)) corresponding 
to Figs. [7] is depicted in Fig. [10] Near the threshold bias, 
the sharp increase of the occupation number of one mode 
is a signature of the runaway solution. We will return to 
the signature in the Raman signal in Sec. |VI| 




FIG. 5: (main panel) Relative magnitude of the BP force 
compared with the average friction 2Fbp/(Ffru +FFR22) as 
a function of level broadening F with ftwo = 0.02 eV. (inset) 
The electronic DOS function at different F. The BP force 
equals the friction at F = 1.2 eV. 
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FIG. 6: Relative magnitude of Berry force compared with 
the average friction felt by the two phonon modes ((Ffru + 
Ffr22)/2) as a function of position of the average electro- 
chemical potential /io at a bias of 1 V. 
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FIG. 7: (Color online) The upper(lower) panels show the 
motion of positive (negative) damped branches calculated at 
V = 2 V bias, with and without the BP force in left and right 
panels, respectively. The middle panel show the inverse Q- 
factor to the left, and the phonon energy as a function of bias 
to the right, with and without the BP force. The solid line 
is the result obtained from the wideband approximation, ig- 
noring the BP force. Parameters used: Fkoo = 20 meV, 8 = 5 
meV. 
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FIG. 8: (Color online) The same as in Fig. [jj but now with 
5 = 4 meV. 
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FIG. 9: (Color online) The threshold bias as a function of 
phonon energy with or without BP force. The parameters 
are the same as in Fig. [jj 
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FIG. 10: (Color online) Semi-classical displacement correla- 
tion function (QQ) of the two modes at different bias. The 
bias increases along the arrow shown, (inset) Excited phonon 
number as a function of bias. The dashed curves are the re- 
sults when only the friction is included. The parameters are 
the same as in Fig. [jj 



is that the ions are moving slowly, while t he ele ctron- 
phonon interaction does not have to be smalP^=. 

In the limit of small ionic velocities, we expand the 
displacement in Eq. ( 10 1 at r' as follows, 



X(t') « Q(t)a z + [ Q(t){t> - t)a z + W)* 



(70) 



where o~ z , I2 are the Pauli and 2x2 identity matrix, re- 
spectively. Using Eq. ( 70 ), we can re-group the expansion 



series in Eq. (10), and arrive at 



G(T,T + )=g (T,T + )+Y] [ g (r,r')M k 



x [Q k (t)(t> - t)a z + ^{t')hj G{T',T + )dr'. (71) 

Now Go(t, t + ) = Qq(t, r+; Q(t)) is the adiabatic electron 
Green's function, determined by the instantaneous elec- 
tronic Hamiltonian when the ions are at a given configu- 
ration (Q(t)). 

The force due to the first term in the new Dyson equa- 
tion now takes the same form as Eqs 



(13 



15| , but the 

non-interacting electron spectral function A a is replaced 
by the adiabatic one, A a (e) = A a (e; Q(t)), which is upto 
infinite order in M, 



Adiabatic limit 



de 
2^ 



TT[A a (e)M k ]An^(e). 



(72) 



The perturbation approach we have presented, and il- 
lustrated with the model calculation above, is applicable 
to weak electron-phonon interaction. It is not restricted 
to slow ions. Within the same theoretical framework and 
based on the Hamiltonian in Eq. (6a|-Eq. (6c), we can 



carry out an adiabatic expansion, where the assumption 



Its contribution to the forces in the Langevin equation 
includes both the renormalization of the effective poten- 
tial and the NC force. To see this, we assume Q is small, 
and expand the adiabatic spectral function over Q near 
Q = 0. The first contribution (Q — 0) is exactly the first 
order result of the perturbation calculation (Eq. (14)). 
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The second contribution (linear in Q) can be split into 
a symmetric and an anti-symmetric part, which are the 
RN and NC force, respectively, 

dQi 

= -2^2 f — ReTr[M k G r {e)M l A a {e)] An a F {e) 
= -2j2 I — ReTr [M k TZ(e)M l A a (e)] An^(e) 
-2J2 I ^lmTi[M k A(e)M l A a (e)] An£(e). 

We can see that in the limit Q — >• they agree with the 
u —> limit of the perturbation results as we should 
expect. Note that above we have treated M as Q- 
independent. If we relieve this assumption, then a further 
term enters the RN forced 

The effective action from the second term in Eq. ( 71 ) 
is given by Eq. (Te)), with Q,(f ) -> Qi{t){t'-t), A a (e) -> 
-4 Q (e). Its contributes to the Langevin equation in terms 
of the friction, the BP force, and the noise. To get the 
expressions for these forces, we need to (1) take the per- 
turbation results in the u —> limit, (2) replace A a (e) 
and Gq(e) with A a {e) and Go (e), respectively. Again, the 
adiabatic results in the Q — > limit agree with the per- 
turbation results in the us — > limit. We conclude this 
section with the diagram showing the relation between 
these two approximations (Fig. 11), and noting that the 



adiabatic approximation in principle allows for updating 
the parameters in the Hamiltonian along the path. 
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FIG. 11: (Color online) The relation between the perturbation 
and the adiabatic expansion. The perturbation expansion 
assumes small deviation from the equilibrium configuration 
(small Q), while the adiabatic approximation assumes slow 
vibrations (small Q). In the region where both Q and Q are 
small, the two expansions agree with each other. 



B. Coupling to electrode phonons 

Actual molecular conductors are coupled also to elec- 
trode phonons. The energy dissipated by the electrons 
can be transferred to the electrodes via this additional 
channeP ^ | 36 | 40 | 79 l. The effective action due to linear 
coupling with a bath of harmonic oscillators is well- 
known 69 If we neglect the electron-phonon interac- 
tion in the electrodes, we can introduce the coupling 
to electrode-phonons in the Langevin equation Eq. ( 28 ) 
by adding the corresponding phonon self-energy: II — 
II e + n p ; l ,II r = W e + 11^. If the phonon baths are at 
equilibrium at a given temperature, they have two effects 
on the system. One is to modify the effective potential, 
and the other is to give rise to dissipation and fluctuating 
forces, which obey the fluctuation-dissipation theorem. 
It is straightforward to include a temperature difference 
between the two phonon baths. A Langevin equation in- 
cluding coupling only with phonon baths has been used 
in molecular dynamics simulations to study phonon heat 



transpor 



It agrees with the Landauer formula in 



the low temperature limit, and with classical molecular 
dynamics in the high temperature limit. 

It is possible to calculate the heat-flux between the 
central region and the phonon reservoirs in the elec- 
trodes using the self-energies describing these baths. If 
we connect the system to the two phonon baths (L and 
i?), and the nonequilibrium electron bath, then the re- 
tarded self-energy giving the deterministic force reads, 
IF = Ilg + W L + U r R , and the fluctuating force is, 
/ = fL + fn + fe- When the system reaches a steady 
state, we can write, 



H p h — I e 



ln = 0, 



(73) 



where I a is the energy current (power) flowing into the 
system from each bath, a (a = L,R,e). Expressions 
for the power exchange can be found using the forces 
acting between the system and each bath in the Langevin 
equation, 

I Q (t) ee -Q T {t) UlT a (t- t')Q(t')dt' - f a (t)j . (74) 

Although we employ the harmonic approximation in this 
paper, this definition is valid also if there is anharmonic 
interaction inside the central regiorP-R This will be 
important for example when describing high-frequency 
molecular modes which only couple via anharmonic in- 
teraction to the low-frequency phonon modes in the elec- 
trodes. This situation can then be handled by a calcu- 
lation where the surface parts of the electrodes are in- 
cluded explicitly in the definition of central region. We 
can write the expression for the energy current in fre- 
quency domain, 



L = (i a (t)) 



(75) 



d/jj 
2^ 



w(tr [U r a {uj)(QQ T ){uj) - (f a Q T )(uj)]) 
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where tr denotes trace over phonon degrees of freedom in 
region C . Using the solution of the Langevin equation, 
and the noise correlation function, we can get a compact 
formula, 



I a = 



duj 
An' 



hu tr 



fl a {uj)D r (u)tl{ui)D a (u) 



Il a (uj)D r {oj)il(uj)D a {uj) 



(76) 



This result agrees with NEGF theorjJal, and fullfills the 
energy conservation, e.g., J2 a =L Re^a =0. Without the 
electron bath and anharmonic couplings, it reduces to the 
Landauer formula for phonon heat transport 83 " 88 !, xhe 
formula contains information about the effects of the elec- 
trons and the electronic current on transport of heat to, 
from and across the central region, but these effects are 
beyond the scope of the present paper. Instead, we now 
focus on how the excitation of the localized vibrations by 
the current affects their Raman signals. 



VI. RAMAN SPECTROSCOPY AND 
CORRELATION FUNCTIONS 

A central aspect of this work is that we need to 
find ways to actually observe the conseque nces of the 
current-induced forces. One promising route^EHl i s t ne 
recent possibility of doing Raman spectroscopy on single, 
current-carrying molecules. In Raman spectroscopy one 
can deduce the "effective temperature" (or the degree 
of excitation) of the various Raman-active vibrational 
modes of a system. The semi-classical theory we have in- 
troduced is not applicable to Raman spectroscopy, since 
it always gives the same Stokes and anti-Stokes lines, 
for a reason that will be clear in the following. So we 
are forced to go back to the quantum-mechanical theory. 
Mathematically, the Raman spectrum can be written as 
follows 

#(w)= [ a k (x k (t)x l (t , ))a l e^ t - t ">d(t~t r ). (77) 



Here a k is a vector involving the change in polarizability 
of the molecule when its atoms are displaced along the 
direction k corresponding to the position operator x k . 
When coupling with the electrode, a k could change due 
to interaction between the molecule and the electrode^!. 
We will take it as a parameter, and focus on the displace- 
ment correlation function instead. Now x k (t) is an opera- 
tor, and the average in Eq. ( |77| is a quantum-mechanical 
one. Since there is no time ordering in the quantum cor- 
relation function (x k (t)xi(t')} at the heart of the Raman 
expression, it is best implemented in our path integral 
version with t' in the upper Keldysh contour and t at the 
lower contour. Hence it can be represented as 



{x k (t)xi{t')) = Z- 



JvQ Jvi \Q k (t) - 



x Qz(f) 



(78) 



where the effective action can be found in formula (16), 

1 f duj 
~2J 2n~ 

-i£t( w )n( w )tf w ) 



Seff(Q,0 



[Qt( w )Lt(«)e(«) + g(u)L(u)Q(u) 

(79) 



and Z is a normalization factor. The Raman spec- 
trum thus has four contributions. A classical con- 
tribution, Mqq(uj), proportional to the average of 
Qh(oS)Qi{uS)* , two quantum corrections, SiQ^(uS) and 
3&£q(uj), proportional to the averages of Qk(ui)£i(u>)* and 
£fc(w)Qi(u;)*, and finally a contribution M^(lo) propor- 
tional to £fc(w)£j(w)*. The calculation of these averages 
involves simple Gaussian integrals, and the results are 

^QqM = a k [hD r {u)tl{uj)D a {uj)] ul ai (80) 



-a k HD r kl (uj)ai 
--a k hD^(uj)ai 



0. 



(81) 

(82) 
(83) 



These functions are dominated by the properties close 
to the poles of L _1 . Let us first consider the case of 
one mode in thermal equilibrium. In this case L(uj) is a 
simple function which can be approximated as 



L(w) 



irjuj. 



(84) 



In the same approximation n(w) is controlled by 
the fluctuation-dissipation theorem and becomes [fi = 

h/{k B T)), 



n(w) = rjui coth 



The classical contribution is now 



up 



T]HiO coth ( ^ 



( — UJ 2 + UJq) 2 + 7] 2 UJ 2 



(85) 



(86) 



This function yields a Raman signal which is symmet- 
ric in uj. It has Lorentzian peaks at u) = ±uj of width 
7]/(2ojq), and with strengths given by its area (integral 
over uj), a 2 Htt / (2ujq) coth(/3a;o/2). Note that the strength 
is proportional to temperature in the high-temperature 
limit. It is also important to note that the strength 
does not depend on the damping, r\. The self-energy, 
fl(w), contributes a factor r\ to the strength, but the oj- 
integration contributes a factor ?7 _1 , hence canceling the 
r] dependence. The physics of this is the fluctuation- 
dissipation theorem of equilibrium: a smaller damping 
should give higher oscillation amplitudes were it not for 
the associated decrease in fluctuations in the environ- 
ment. 

The quantum correction ^q^(lo) + &£q(w) is propor- 
tional to the imaginary part of L _1 . In the above ap- 
proximation it becomes, 



M Qii {u) + M iQ {u) = a 2 



rjhui 



{—uj 2 + ojq) 2 + rfuj 2 



(87) 
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This contribution breaks the uj — u symmetry, and 
is hence responsible for the different strengths of the 
Stokes(phonon emission) and anti-Stokes(phonon ab- 
sorption) lines. This term also exhibits peaks at ±a>o, 
with strengths ±a 2 Htt / (2wo) . 

The ratio, r, of strengths of the anti-Stokes and Stokes 
lines becomes, 



coth(^) 
cothte 



In equilibrium this ratio is used to measure the temper- 
ature employing the particular mode excitation, but has 
also been used to estimate an effect ive te mperature of 
the different modes out of equilibriu m 1 89 * 93 !. 

Next we discuss the situation in the presence of elec- 
tronic current and the derived forces. We consider a sit- 
uation with two modes close in frequency, which are cou- 
pled by the current. Now, the function L(uj) is a 2 x 2 
matrix. If the NC and BP forces are represented by con- 
stants, n and b, respectively, then L(lo) becomes 



L(w) 



+ U>1 — 11x17] 

-n + iuib 



n — iujb 



(89) 



Here we have made the simplifying assumption that the 
friction, r), is the same for the two modes. The two mode 
frequencies uj\ and uji have a difference 5 = u>2 — Wi- 



Comparing with the expressions Eqs. (57 1, (59), the pa- 



rameters n and b are both linearly dependent on the ap- 
plied voltage V. If this voltage is sufficiently large, there 
is a possibility that one of the eigenmodes of L(ui) will 
have a vanishing imaginary part at a critical voltage V c , 
and we obtain a runaway mode. 

In the following we present plots of the resulting Ra- 
man spectra as a function of voltage approaching V c 
from below. The function II(w) describing the fluctuat- 
ing forces is both temperature- and voltage-dependent. 
However, nothing dramatic happens at the critical volt- 
age, so it will be taken to be a temperature and voltage 
independent matrix with values of the order of rju. As in 
the equilibrium case, an increasing lifetime of the mode 
will lead to stronger oscillation amplitudes, but out of 
equilibrium this need not be counteracted by decreasing 
fluctuations in the environment. The result is a strong 
increase in the strength of the runaway mode, both for 
the Stokes and the anti-Stokes line. 

Figure [12] shows the Stokes lines for two modes close in 
frequency and coupled by the NC and BP forces, for the 



system described by Eq. ( 89 ) . In the following graphs the 
parameters are, lo\ — uj—5/2, u>2 — uj+S/2 with 5 = 0.2w, 
7/ = 0.05w, b/n = 0.6/oj, for which the critical value 
of the strength of the NC force will be n c = 0.124cj 2 , 
when the voltage reaches its critical the value, V c . We 
see that the mode at L02 is increasing in strength while 
its frequency is shifting slightly down. By contrast the 
frequency of the mode at oji is moving upwards while the 
strength is decreasing. If we compare to the standard 



one-mode, equilibrium theory of Raman lines, we would 
say that one mode heats up, while the other cools down. 
This is seen clearly in Fig. [l3j where the strength of the 
two modes is plotted as a function of voltage. Here we 
see a small cooling of one mode, while the other mode 
has a diverging temperature, as the voltage approaches 
V c . Alternatively, one could use the ratio of the anti- 



Stokes/Stokes strengths in combination with Eq. (88) to 



determine an effective mode-temperature, as shown in 
Fig. [14] Interestingly, we cannot see the cooling effect 
from the anti-Stokes/Stokes ratio. 

This could be qualitatively understood as follows. In 
principle, Eq. (88) only holds at equilibrium, which is 
a result of the fluctuation-dissipation theorem. Under 
nonequilibrium conditions, discrepancies between differ- 
ent ways of defining an effective temperature may be ex- 
pected. Specifically, in the case studied here, the bias 
modifies the rates for phonon emission and absorption. 
This can be inferred from the change of peak broadening 
in the phonon correlation function, as well as from the 
shift of the peaks, with bias. Part of this bias-dependent 
effect is lost, when we form the anti-Stokes/Stokes ratio. 
But it is included when we look at the bias-dependent 
strength of the Stokes lines. 




0J\ (jj 2 

Raman Frequency 



FIG. 12: (Color online) The Stokes lines of two coupled modes 
for various values of the applied voltage in units of the critical 
voltage. Parameter values are u)\ = uj — 8/2, u>2 — oj + 8/2 
with 8 = 0.2lj, 77 = 0.05o>, b/n = 0.6/lj. 



VII. CONCLUDING REMARKS 

In this paper we have derived a semi-classical Langevin 
equation describing the motion of the ions, in the har- 
monic approximation, in nanoscale conductors including 
both the effective action from the current-carrying elec- 
trons and the coupling to the phonon baths in the elec- 
trodes. Joule heating and current-induced forces are de- 
scribed on an equal footing by this methodology. We 
derive a convergent expression for the BP force, remov- 
ing the infrared divergence in our previous result!^. The 
importance of the BP force in relation to the stabil- 
ity of the device is further highlighted in a two-level 
model system. Using the same model we show the signa- 
ture of the current-induced runaway mode excitation in 
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FIG. 13: (Color online) Strength of the two Stokes lines as a 
function of voltage. Solid line corresponds to the "runaway" 
mode at L02, and dashed line to the "cooling" mode at Wi. 




FIG. 14: (Color online) Ratio of strengths of the anti-Stokes 
and Stokes lines as a function of voltage. Parameters are as 
in Fig. [12] Solid line corresponds to the "runaway" mode at 
u)2, and dashed line to the "cooling" mode at wi. 



the Ram an spe ctroscopy of a current-carrying molecular 
conducto r 89 * 90 !. 

We should mention that the harmonic approximation 
used here breaks down when the phonons are highly ex- 
cited. In that case the anharmonic coupling between dif- 
ferent phonon modes becomes important. Analysis in 
this regime relies e.g. on molecular dynamics simulation. 
To take the nonlinear potential approximately into ac- 
count in the present method, we only need to replace the 
force from the harmonic potential —KQ(t) with the full 
potential — 8qV(Q). Thus the full Langevin equation has 
the following advantages: first, we can include the highly 
nonlinear ion potential, which is necessary to simulate 
bond-breaking processes; second, a crucial part of the 
quantum-mechanical motion is included in the dynamics. 
For example, we recover the quantum-mechanical results 
for Joule heating and for heat transport in the harmonic 
limit. This enables the study of phonon heat transport!^, 
or thermoelectric^ transport using the same Langevin 
equation including the electron-phonon interaction, and 
is an interesting topic for future research. 



The backaction of the runaway vibrations on the elec- 
trons is another possible extension of the present study. 
There are at least two effects to address. The first one 
is the adiabatic change of the electronic structure. The 
purpose of the adiabatic extension in Sec. V A is to in- 
clude this effect. The second effect is the inelastic electri- 
cal current, which becomes important when the electron 
mean free path is comparable to the device length. 

A number of problems need to be faced when imple- 
menting the approach within the framework of density 
functional theory. First, the ions may be driven away 
from their equilibrium positions by the current. The elec- 
tronic structure and electron-phonon coupling depend on 
the ionic positions, and one may need to update the elec- 
tronic friction and noise correlation function throughout 
the molecular dynamics simulation. This introduces a 
technical problem, in addition to the computational chal- 
lenge, which is how to generate colored noise when its 
correlation function is time-dependent. Second, calcula- 
tion of the convolution kernel in the Langevin equation is 
time-consuming considering the many time-steps needed 
to sample the dynamics. For the electron bath in the 
wideband limit the convolution transforms to time-local 
forces. But the time-scale of the phonon bath is typi- 
cally comparable to that of the system, and we cannot 
use the wideband approximation. One possible solution 
could be to include more ions from the baths into the 
dynamical region and approximate the coupling to the 
external phonon bath with time-local forces. Intuitively, 
with more ions included, the approximate central system 
will be closer to the actual system under study. Work to 
overcome these difficulties is underway. 

Electron-nuclear dynamics is a broad problem, relevant 
to many fields. A central challenge is how to take ac- 
count of electron-nuclear correlation. The simplest form 
of nonadiabatic dynamics - the Ehrenfest approximation 
- fails precisely there: it does not take into account spon- 
taneous phonon emission by excited electrons and conse- 
quently the Joule heating effect. The appeal of Ehren- 
fest dynamics is its conceptual simplicity derived from 
the classical treatment of the nuclei. The Langevin ap- 
proach retains this key element, while rigorously reinstat- 
ing the vital missing ingredient: the fluctuating forces - 
with the correct noise spectrum - exerted by the electron 
gas and responsible for the return of energy from excited 
electrons to thermal vibrations. Thus, in addition to its 
capabilities as a method for nonadiabatic dynamics, this 
approach can be very helpful conceptually, by explicitly 
quantifying effects that are intuitive but whose physical 
content can sometimes remain hidden from view. 

Although our discussion here is in the context of 
molecular electronics, we believe that the predictions 
based on the simple model can be important for 
other interesting physical systems, for instance nano- 
electromechanical o scillat ors (NEMS) coupling with an 
atomic point contadP^l When the current through an 
atomic point contact is used to measure the motion of 
an oscillator, the measurement imposes quantum back- 
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action on the oscillator. A similar H amilton ian describes 
this quantum backaction quite wel l 61 * 67 * 98 *. If we now 
couple two identical oscillators, or two nearly-degenerate 
eigenmodes by the point contact, it seems to be experi- 
me ntally f easible to detect the polarized motion predicted 
herepEM. 

A neighbouring field, where much larger size- and time- 
scales come into play, and where Langevin dynamics is an 
import ant line of approach, is the simulation of radiation 
damage^M It is hoped that the present discussion will 
be of interest to that community. 



In the wideband limit, ignoring the w dependence of 
ImAA, we get 



B 



2heV 

7tO, 



-x 



(A5) 



where fl c is an upper bound on the electron-hole pair ex- 
citation. In Ref.l55( we used the largest phonon frequency 
instead, which over-estimates the effect of BP force. If as 
a conservative estimate of f2 c we take a typical hopping 
matrix element, then using Q/Q ~wwe get the estimate 
in ([541). 
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Appendix A: Connection with Nano Lett., 10, 1657 
(2010) 

In Ref. 1551 we carried out an adiabatic expansion in a 



slightly different way from that in Sec. V A We obtained 
an infrared divergence in the expression for the BP force, 
and used the largest phonon frequency as a cutoff. Here 
we show that by introducing a small lifetime broaden- 
ing to the scattering eigenstate (7), we can remove the 
divergence, and get the same result as shown in this pa- 
per. We start from the expression for the Berry force in 
Ref. [55J an d write it as 

B = - hm / dw- (Al) 

7->oy (uj-ij) 2 

We first do a partial integration to get 



B 



lim 

7->0 J U> — 17 



^ImAAfa;) 



(A2) 



From the w-dependent part of ImAA(w)(Eq. (45 1) 



7->0 J uj — vy 

(A3) 

which gives 

B pa 2h^2 f — An£(e) 

x lmTi[MA a (e)Md e Ren(E)}. (A4) 
This agrees with the result in Sec. |V A 



Appendix B: Alternative derivation of the Langevin 
equation 

In this Appendix, we give an alternative way of arriv- 



ing at the generalized Langevin equation (Eq. (28)). The 



derivation here is meant to be intuitive rather than the- 
oretically rigorous. Our starting point is the equation of 
motion for the (mass-normalised) displacement operator 
u. 



x = —Kx + F e , 
where we define the electronic force operator 

dH e (x) 



F, 



dx 



(Bl) 



(B2) 



With the help of Green's functions, Eq. (Bl I can be cast 



into the following form, for each degree of freedom k, 

x k = - KkjXj + ihTr [M k G < {t, t+)] + f k (tXB3) 
3 

We have defined the noise operator 

h{t) = iJ2 M L (<4(*M*) - ftG< n (M+)) , (B4) 



and the lesser Green's function is G< m (t,t + ) — 
(«//i)(cj n (t + )c n (t)). The quantum average (...) is over 
the electronic environment, which need not be in equilib- 
rium. 



From Eq. (10) to second order in M, 



G<(t,t + ) = G<(t,t+) 

[ G (t,t')M k x k (t')G<(t',t + )dt' 
k J 

J G<{t,t')M k x k (t')G {t\t+)dt' . 



(B5) 



Using this in Eq. (|B3|), we get 



kj x j 



iHTr[M lc G$(t,t + )} 



J Tl r kj (t,t')x 3 (t')dt' + f k (t), (B6) 
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which is of the same form as Eq. ( 28 1 



Now consider the time correlation of the noise operator 



and 



(fi(t)fj(t')) = h 2 TT[M l G>(t,t')M^G<(t',t)} 
= ^MI> (<,*'), 



(f J (t , )Mt))=thu<(t,t / ). 



(B7) 



(B8) 



As expected, the quantum-mechanical noise operators at 
different times does not commute. 

To go to the semi-classical approximation, we take the 
classical noise correlation as the average of Eqs. (B7) and 
(PI 



m)f j (t , ))c = (f j (ty i (t)) c 



(B9) 



where (. . . ) c denotes a classical statistical average. Now 
the noise spectrum becomes "classical" , in the sense that 
its time correlation function is real. If the potential is 
anharmonic, the right side of Eq. (B6 1 will contain terms 



of higher order in u. The equations of motion of these 
higher-order terms form an infinite hierarchy. Then the 
semi-classical approximation will have to involve a trun- 
cation procedure. However, this is a separate problem, 
beyond the scope of this paper. 

We stress again that after making the semi-classical ap- 
proximation, it is still possible to calculate the quantum- 
mechanical average of two displacement operators at 
equal times (tii(t)uj(t)) within the harmonic approxi- 
mation, from the semi-classical Langevin equation, lead- 
ing to the correct quantum-mechanical vibrational energy 
and steady-state transport properties. 



Appendix C: The phonon self-energy 

The phonon self-energies correspond to the bubble di- 



agram in Fig. [T] 



n 



/,-/ 



(t-t') = -ih r Tt[M k G^' > {t-t')M l G^ < {t' -t)], 
nfci°(*-0 = -ifrTr[M k G r Q ' a (t-t')M l G£(t' -t)] 

-iftTr [M k G<(t - t')M l G% r (t' - t)](.Cl) 
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